#This calculates temperature and resistance using the ITS-90 calibration standard 
# for a platinum RTD. 
#For reference see www.its-90.com
from math import exp,log,fsum,pow

#ITS-90 constants
WR={}
T={}
WR['H2O']=1.00000000
T['H2O']=0.01000
WR['Ga']=1.11813889
T['Ga']=29.7646
WR['Hg']=0.84414211
T['Hg']=-38.8344
WR['Ar']=0.21585975
T['Ar']=-189.3442
WR['In']=1.60980185
T['In']=156.5985
WR['Sn']=1.89279768
T['Sn']=231.928

class calibration(object):
	def __init__(self):
		#Standard set of parameters, perfect calibration
		self.set(low=[0,0],high=[0,0],R0=100,T0=273.16,Roffset=0.0)
		
	def set(self,**dicty):
		for key,value in dicty.items():
			if (key not in dir(self)) and (key.lower() in [val.lower() for val in dir(self)]):
				key=next(val for val in dir(self) if val.lower()==key.lower())
			self.__setattr__(key,value)
		
	@property
	def R0(self):
		return self._R0
	@R0.setter
	def R0(self,value):
		self._R0=value

	@property
	def low(self):
		return self._low
	@low.setter
	def low(self,value):
		self._low=list(value)
	
	@property
	def high(self):
		return self._high
	@high.setter
	def high(self,value):
		value=list(value)
		if len(value)<3: value.extend([0.0,0.0,0.0])
		self._high=value
		
	@property
	def a_low(self):
		return self.low[0]
	@a_low.setter
	def a_low(self,value):
		self.low[0]=value
	@property
	def a_high(self):
		return self.high[0]
	@a_high.setter
	def a_high(self,value):
		self.high[0]=value
	@property
	def b_low(self):
		return self.low[1]
	@b_low.setter
	def b_low(self,value):
		self.low[1]=value
	@property
	def b_high(self):
		return self.high[1]
	@b_high.setter
	def b_high(self,value):
		self.high[1]=value
		
	@property
	def Range(self):
		# returns a list of min,max temperature in units of K
		if ('_range' not in dir(self)):
			#defaults to melting pt of mercury as low, melting pt of Gallium as high
			return [234.3156,302.9146]
		return self._range
	@Range.setter
	def Range(self,value):
		# requires input in units of K
		if ('_range' not in dir(self)): 
			self._range=[]
		if type(value) is float: self._range.append(value)
		if type(value) is list: self._range=value

class RTD(object):
	R_units='ohm'
	T_units='K'
	
	A=(-2.13534729,3.18324720,-1.80143597,0.71727204,0.50344027,-0.61899395,-0.05332322,\
		0.28021362,0.10715224,-0.29302865,0.04459872,0.11868632,-0.05248134)
	enumA=tuple(enumerate(A))
	
	B=(0.183324722,0.240975303,0.209108771,0.190439972,0.142648498,0.077993465,0.012475611, \
		-0.032267127,-0.075291522,-0.056470670,0.076201285,0.123893204,-0.029201193,-0.091173542,\
		0.001317696,0.026025526)
	enumB=tuple(enumerate(B))
		
	C=(2.78157254,1.64650916,-0.13714390,-0.00649767,-0.00234444,0.00511868,0.00187982,\
		-0.00204472,-0.00046122,0.00045724)
	enumC=tuple(enumerate(C))
		
	D=(439.932854,472.418020,37.684494,7.472018,2.920828,0.005184,-0.963864,-0.188732,\
		0.191203,0.049025)
	enumD=tuple(enumerate(D))
	del(A,B,C,D)
	
	def __init__(self):
		self.Tout=True
		self.units='C'
		self.cal=calibration()
		
	def __call__(self,R):
		self.R=R
		return self.T
	
	@property
	def units(self):
		return self._units
	@units.setter
	def units(self,value):
		self._units=value
		K=lambda x: x
		C=lambda x: x-273.15
		F=lambda x: C(x)*9./5.+32
		R=lambda x: x*9./5.
		if self.units=='K': self.convert=K
		if self.units=='C': self.convert=C
		if self.units=='F': self.convert=F
		if self.units=='R': self.convert=R
		
	@property
	def R(self):
		#returns in units of Ohms
		if self.Tout:
			return self._resistance
		return self.W*self.cal.R0
	@R.setter
	def R(self,value):
		#requires input in units of Ohms
		self.Tout=True
		self._resistance=value+self.cal.Roffset
	@R.deleter
	def R(self):
		self.Tout=False
		
	@property
	def T(self):
		#returns in units of as defined in self.units
		if self.Tout:
			if (self.R>self.cal.R0):
				return self.convert(self.eq10b())
			else:
				return self.convert(self.eq9b())
		return self._temperature
	@T.setter
	def T(self,value):
		#requires input in units of K
		self.Tout=False
		self._temperature=value
	@T.deleter
	def T(self):
		self.Tout=True
	
	@property
	def W(self):
		# W is the ratio of measured resistance to resistance at the triple point of water. 
		# It is unitless.
		if self.Tout:
			return self.R/self.cal.R0
		else:
			#Add b instances of 13 and 14 that solve for W from Wr
			#Then reference here depending on temperature
			pass
			
	@property
	def Wr(self):
		# Wr is the deviation function. It is unitless
		if self.Tout:
			if (self.R>self.cal.R0):
				return self.eq14a()
			else:
				return self.eq13a()
		else:
			if (self.T>self.cal.T0):
				return self.eq10a()
			else:
				return self.eq9a()
		return None
		
	def eq13a(self):
		# Valid over the range 83.3K < T90 < 273.16K
		# This function returns the deviation function. It is unitless
		W=self.W
		C=(self.cal.low)
		return fsum([-C[0]*(W-1),-C[1]*(W-1)*log(W),W]) #Wr
				
	def eq14a(self):
		# Valid over the range 0C < T90 < 961.78C
		# This function returns the deviation function. It is unitless
		W=self.W
		value=W-1.0
		C=self.cal.high
		return fsum([-C[0]*value,(-C[1]*value**2.),(-C[2]*value**3.),W]) #Wr

	def eq9a(self):
		# This function valid over the range 13.3K < T90 < 273.16K
		# This function returns the deviation function
		value=((log(self.T/273.16)+1.5)/1.5)
		return exp(fsum([Coef*pow(value,i) for (i,Coef) in self.enumA])) #Wr
		
	def eq9b(self):
		# This function valid over the range 13.3K < T90 < 273.16K
		# This function returns temperature in K
		value=((self.Wr**(1/6)-0.65)/0.35)
		return fsum([Coef*pow(value,i) for (i,Coef) in self.enumB])*273.16 #T
		
	def eq10a(self):
		# This function valid over the range 0C < T90 < 961.78C
		# This function returns the deviation function
		value=((self.T-754.15)/481)
		return fsum([Coef*pow(value,i) for (i,Coef) in self.enumC]) #Wr
		
	def eq10b(self):
		# This function valid over the range 0C < T90 < 961.78C
		# This function returns temperature in K
		value=((self.Wr-2.64)/1.64)
		return fsum([Coef*pow(value,i) for (i,Coef) in self.enumD])+273.15 #T

def calibrate(R0, **points):
	from numpy.linalg import solve
	points=[calibration_point(point) for point in points.values()]
	RTD_Ws=[]
	RTD_Wrs=[]
	item=RTD()
	item.cal.R0=R0
	for point in points:
		del(item.R)
		item.T=point.T
		RTD_Wrs.append(item.Wr)
		del(item.T)
		item.R=point.R
		RTD_Ws.append(item.W)
	#print(RTD_Ws, RTD_Wrs)
	if all(point.T>273.15 for point in points):
		A=[[(RTD_Ws[0]-1),((RTD_Ws[0]-1)**2)],[(RTD_Ws[1]-1),((RTD_Ws[1]-1)**2)],]
	else:
		A=[[(RTD_Ws[0]-1),((RTD_Ws[0]-1)*log(RTD_Ws[0]))],[(RTD_Ws[1]-1),((RTD_Ws[1]-1)*log(RTD_Ws[1]))],]
	B=[(RTD_Ws[0]-RTD_Wrs[0]),(RTD_Ws[1]-RTD_Wrs[1]),]
	return solve(A,B)

class calibration_point(tuple):
	@property
	def T(self):
		return self[0]
	@property
	def R(self):
		return self[1]
		
		
		

import unittest
class Perfecttest(unittest.TestCase):
	def setUp(self):
		self.RTD=RTD()
	def keycheck(self,key):
		self.assertAlmostEqual(self.RTD(WR[key]*100),T[key],places=3)
		self.assertAlmostEqual(self.RTD.Wr,WR[key])
		self.assertAlmostEqual(self.RTD.W,WR[key])
	def test_Ar(self):
		self.keycheck('Ar')
	def test_Hg(self):
		self.keycheck('Hg')
	def test_H20(self):
		self.keycheck('H2O')
	def test_Ga(self):
		self.keycheck('Ga')
	def test_In(self):
		self.keycheck('In')
	def test_Sn(self):
		self.keycheck('Sn')

if __name__=='__main__':
	instr=RTD()
	print('test out calibrations: M240160')
	M240160=RTD()
	calpoints={'R0':99.98}
	calpoints['point1']=(39.98635+273.15,115.47368)
	calpoints['point2']=(99.94699+273.15,138.42492)
	M240160.cal.high=(calibrate(**calpoints))
	print('high calibration coefficients:',calibrate(**calpoints))
	calpoints['point1']=(-9.98543+273.15,96.13693)
	calpoints['point2']=(-19.98702+273.15,92.22590)
	M240160.cal.low=(calibrate(**calpoints))
	print('low calibration coefficients: ',calibrate(**calpoints))
	M240160.cal.set(R0=calpoints['R0'])
	
	res=[189.279768, 138.42492,126.970076,115.47368,100.00,99.97833,96.13693,92.22590,21.585975]
	print('Calibrated temperatures are:')
	[print('R={:8.4f}, T={:8.4f}'.format(value,M240160(value))) for value in res]
	import timeit
	test1=timeit.Timer('[M240160(value) for value in res]',"from __main__ import M240160, res")
	n=10000
	print('timing for {} calls in the following cases:'.format(n))
	print('object call timing for temperature:    ',test1.timeit(n))
	
	

	
	